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Abstract. Recent experiments have shown that various structures may be formed 
during the evaporative dewetting of thin films of colloidal suspensions. Nano-particle 
deposits of strongly branched 'flower-like', labyrinthine and network structures are 
observed. They are caused by the diff'erent transport processes and the rich phase 
behaviour of the system. We develop a model for the system, based on a dynamical 
density functional theory, which reproduces these structures. The model is employed to 
determine the influences of the solvent evaporation and of the diffusion of the colloidal 
particles and of the liquid over the surface. Finally, we investigate the conditions 
needed for 'liquid-particle' phase separation to occur and discuss its effect on the self- 
organised nano-structures. 



Modelling the evaporation of thin films of colloidal suspensions using DDFT 2 
1. Introduction 

Surface patterns resulting from structure formation occur naturally in many different 
systems and are extensively studied in various scientific fields. A classic example are 
branched patterns, e.g. found in river networks [Ij or formed by bacterial colonies 
[2], that sometimes form fractals. Other examples are the labyrinth patterns such as 
those formed via calcification and mineralisation processes [3j. The spontaneous self- 
assembly and self-organisation of atoms, molecules and nano-particles at interfaces is 
a widely researched topic not only because the resulting structures are interesting but 
also because they may be used in the manufacturing of nanostructures, e.g. assembling 
colloids to form photonic bandgap crystals [4j. A large variety of intricate structures 
can be formed even during the dewetting process of films of non- volatile fiuids on solid 
substrates. The dewetting process starts with the rupture of the initially homogeneous 
film caused either by a surface instability (spinodal dewetting [5j) or by the nucleation 
of holes which often occurs at surface defects [HI [71 [HI |9l [JOj [HI [I2j . The holes then grow 
[13] and the rims meet to form a global network, drop or labyrinth pattern p!4[ [T5l [16]. 

The evaporative dewetting of polymer/macromolecule solutions [TTl HU [I9] and 
colloidal (nanoparticle) suspensions [20l [211 (ZSl [23l [24] can produce a wide variety of 
patterns and has been intensely studied in various experimental settings over recent 
years. Although the motivation for the work presented here is mainly drawn from 
the latter, we believe that our results also explain the basic features in the case of 
evaporative dewetting of solutions. The particular experiments that directly motivate 
our theoretical work presented here are those described in Refs. [I22[ [24l [25l |26l [23] that 
use a suspension of thiol-passivated gold nano-particles dispersed in an organic solvent. 
A drop of the suspension is spin-coated onto a fiat silicon substrate to form a thin film 
over the surface. The solvent evaporates during the spin-coating and leaves a nano- 
particle pattern on the surface. In another experimental setup a drop of suspension 
is placed on the surface within a tefion ring [27j. The evaporative dewetting is slower 
than in the case of spin-coating and the structuring is observed using video-microscopy 
p3j . What these experiments show is that branched structures are formed by transverse 
instabilities of the receding mesoscopic contact line. However, more intriguing are the 
patterns that are formed in an ultrathin layer that is left behind this contact line. 
Three different types of structures have been observed: a labyrinth pattern formed 
during spinodal dewetting, a two-scaled network structure formed via the nucleation 
and growth of holes and a branched structure formed by a fingering instability that 
occurs at the dewetting front of nucleated holes. The diffusive mobility of the nano- 
particles and the interaction energies between the particles can be altered by changing 
the thiol chain length. The fingering instability is found for relatively low nano-particle 
mobilities. Because these structures form in the ultra-thin layer, the height of film is of 
the same order of magnitude as the diameter of the colloids. 

To model these processes, simple two-dimensional (2d) lattice kinetic Monte-Carlo 
(KMC) models were proposed j28l |26l [23l [^J. Vancea et. al. [29] made a detailed 
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investigation of the characteristics of the branched structures and their dependence 
on the model parameters. They observed that as in the experiments the fingering 
instabihty becomes stronger the smaUer the nano-particle mobihty gets. A pseudo 
three-dimensional KMC model has also been considered [261 [30 ] which can reproduce 
dual-scale network patterns. 

An alternative approach that may be used for modelling such systems is based 
on thin-film hydrodynamical models, which are derived by making a long-wave 
approximation [31j. Recently, for example, line-pattern formation has been observed in 
a simple long-wave model for a thin film of colloidal suspension evaporatively dewetting 
from a surface [32j. Such thin film equation based models provide a good description 
of the system on mesoscopic length scales and reproduce the experimental results 
[331 [341 l35l [36j. However, they are unable to describe the dynamics of the system 
at the microscopic (single particle) level. In Ref. [37j a more detailed account of 
the different approaches that may be used for modelling the evaporative dewetting of 
colloidal suspensions is given, so we do not re-review the subject here. 

In the present work, we give a detailed account of an alternative model for this 
system that has been briefiy discussed before [371 EH] . It is based on Dynamical Density 
Functional Theory (DDFT) [39l|40l [HI |42] and goes beyond the 2d KMC model by 
also allowing us to investigate the infiuence of liquid diffusion over the surface. This 
paper is laid out as follows: First, in Sec. [2} we present the coarse-grained model for the 
Hamiltonian and the resulting approximation for the free energy used in our model. This 
is followed in Sec.|3]by an introduction to the dynamical equations of the DDFT. In Sec. [4] 
we discuss the phase diagram and perform a linear stability analysis of homogeneous 
steady films, whereas in Sec. [5] we present fully nonlinear simulation results. In Sec. [6] 
we summarise our findings and draw some conclusions. 

2. Free Energy for the System 

We consider a simplified coarse-grained two-dimensional model for the system. The 
surface of the substrate is divided into a square array of discrete lattice sites. We 
choose the cell size so that each cell may be occupied by at most one nano-particle, 
i.e. the lattice spacing a is roughly equal to the diameter of the nano-particles. As 
a consequence, each lattice site must be in one of three possible states: (i) occupied 
by a nano-particle, (ii) occupied by liquid or (iii) unoccupied (as displayed in Fig. [T]). 
To characterise the state of the system, we introduce two occupation numbers for each 
lattice site i\ Ui for nano-particles and U for liquid. These occupation numbers are 
binary values which describe the state of each site. The occupation numbers for (i) a 
lattice site containing a nano-particle would be = 1, = 0, for (ii) a site occupied by 
a film of liquid we have = 0, = 1 and for (iii) a vacant site we have = 0, = 0. 
This type of lattice model has the following Hamiltonian [28]: 



E = — ei ^ lilj — Cn 




<ij> <ij> 



<ij> 




Figure 1. The sketch indicates how the substrate is divided into lattice sites and 
shows the three possible states of each lattice site. 



i i i 

where J2<ij> denotes a sum over pairs of nearest neighbours. (/)[ and 0^ are external 
potentials acting on the liquid and nano-particles respectively, at lattice site i. Three 
interaction terms are included which determine the strength of attraction between 
neighbouring cells, ei is the interaction energy between two adjacent cells containing 
films of liquid, is for adjacent cells which both contain nano-particles and Cni is 
the energy between a cell containing a nano-particle and a cell containing liquid. The 
amount of liquid on the surface is not a conserved quantity because it can evaporate to 
and condense from a reservoir of vapour above the surface, /x is the chemical potential 
of this reservoir. 

From the Hamiltonian ([T]) we can derive an expression for the free energy of the 
system. The probability that a lattice site i is covered by a liquid film in an equilibrium 
configuration is given by the following integral over time t: 

1 r 

p\ = lim - / k{t)dt. (2) 

Similarly, the probability that a lattice site i contains a nano-particle is given by the 
following expression: 

= lim - / ni{t)dt. (3) 

By choosing the grid spacing a to be our unit of length, equal to one, these probabilities 
are also equal to the local number densities for the liquid and nano-particles. Following 
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the approach described in Refs. |43l [44j . we make a Bragg- WiUiams mean field 
approximation for the (semi-grand) free energy of the system: 

F = kBTY, [P\ In + (1 - - p^) + pr In p^ 



+ (l-pnin(l-pn]-e;E^M 



<ij> 



(4) 



where ks is Boltzmann's constant and T is the temperature. This is a semi-grand 
free energy because the hquid is treated in the grand canonical ensemble (the reservoir 
of vapour fixes the chemical potential /x), whereas the number of nano-particles in the 
system is a conserved quantity (so these are treated canonically) . Eq. Q is derived using 
a 'zeroth-order' mean-field approximation and thus higher order terms are omitted (e.g. 
the terms involving ln(l — — p^) which describe the excluded area correlations between 
the nano-particles and the liquid). If we assume that the density values only vary on 
length scales 3> a, then we can define a gradient operator for this discrete system: 

y+l - Pl,y) ■> 

^ Px,y — iPx+l,y ~ Px,yi Px,y+1 ~ Px,y ), (5) 

where each lattice site i is now represented by the coordinates (x^y). Using the operators 
([5]) we can express the summation over pairs of nearest neighbours as 

E p?p' = 4 E f^p' - Y.(^p?) ■ (vpf ), (6) 



<ij> 



where f3 = I. Substituting ^ into our lattice free energy Q we obtain 
F = kBTj2[p\ Inp^ + (1 - p^) ln(l - p^) + p^ Inp^ 

i 

+ (1 - pn in(i - pd] - E \U(p^' + t'-^p'^f + ^'-ipy\ 



+ 



^(Vp[f + f(Vprf + e„,(Vpn-(Vpi) 



(7) 



where the factor of \ is included to avoid double counting. Taking the continuum limit, 



so 



that Jd^^ p7 Pn{r), p\ Pi{y), 



(r) and 



(r) where the 



vector r = (x, y) is a continuous variable, we obtain for the free energy of the system: 



dr 



/(p,(r),p.(r)) + |(VpKr))2 



+ |(Vp„(r)f + e„,(Vp„(r)).(Vp,(r)) 



+ 



J drp;(r)(^i(r) - p) + J (/rp„(r)^„(r), 



(8) 
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where 



f{ph Pn) = kBT[pi \jipi + (1 - pi) ln(l - pi) + Pn \VLPn 

+ (1 - Pn) ln(l - Pn)] - 2eip^ - 2enpl " 4e^/PnA. (9) 
This free energy functional may be employed to determine the phase diagram of the 
system - i.e. the state of the system in the thermodynamic limit (see Section [4]). 
However, the observed patterns are often non-equilibrium structures that are 'dried in', 
i.e., that evolve towards the equilibrium state on a time scale that is much longer than 
the typical observation times. To model the non-equilibrium processes that result in the 
observed self-organised structures, one needs kinetic equations for the time evolution of 
the densities. They are developed in the following section. 



3. Modelling the Dynamics of the System 



The chemical potential of the nano-particles may be calculated using the following 
functional derivative 051 [46l [47] : 

SF[pn,pl] 



Sprj 



(10) 



In equilibrium systems the chemical potentials take a uniform value throughout the 
system. However, this is not the case for non-equilibrium configurations that the system 
takes during its time evolution. There, the chemical potential varies temporally and 
spatially over the surface. In particular, non-equilibrium density profiles pi{r^t) and 



Pn{v^t) give, via Eq. (10), a non-equilibrium chemical potential for the nano-particles 
/i^(r,t). Thus, the time-dependent densities are coarse-grained 'average' quantities. 
Assuming that locally, the system is in equilibrium, we may define these non-equilibrium 
density fields in a similar way as the equilibrium densities [see Eqs. ^ and ([3])]: 

p\ = — km. (11) 

tm Jo 

p1 = — ni{t)dt, (12) 
tm Jo 

where tm is now a finite time that is large compared to the solvent molecular collision 
time, but is small compared to the time scale for a nano-particle to move from one 
lattice site to a neighbouring lattice site. If we now assume that the driving force which 
causes a flux of the nano-particles over the surface of the substrate to be given by the 
gradient of the chemical potential then the nano-particle current is given by 



-M,(r,t)V- 



(13) 



where M„(r, t) is a mobility coefficient which we assume to depend on the local densities 
Pn{r,t) and pi(r,t). Since the number of nano-particles in the system is conserved we 



can combine Eq. (13) with the continuity equation to get 

,5F[pn,piy 



dpn{r,t) 
dt 



Mn{Pn,Pl)^- 



Spnir,t) 



(14) 
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We expect that when the hquid density is uniform throughout the system {pi{r) = p) 
and the density of the nano-particles is smaU everywhere <C 1) then the nano- 



particle dynamics given by Eq. (14) must reduce to the diffusion equation (i.e. Eq. (13) 



becomes Pick's law). When we apply this to Eq. (10), for small p^, we get the leading 
order term: 



/i,(r) ^ kBT\liPn{v) + C + 0{pn), 



(15) 



where C = C{pi) represents constant terms. Substituting Eq. (15) into Eq. (14) we see 



that in order for Eq. (14) to reduce to the diffusion equation the mobility has to be 



proportional to Pn{v^t). Therefore, we write the nano-particle mobility as: 

Mn = Pn{v,t)m{pi{v,t)), (16) 

where m{pi{r^ t)) is a function of the local liquid density. We assume that the diffusion of 
the nano-particles over the surface is caused by the Brownian 'kicks' from the molecules 
in the liquid. Therefore when the liquid density is low (on the dry substrate) the nano- 
particles are almost immobile. However, when the density of the liquid on the surface 
is high (where the substrate is covered by the liquid film), the nano-particles are much 
more mobile. We model this behaviour by a function m(p^(r, t)) which switches between 
a very small value, when pi is small and a, which is the mobility coefficient for the nano- 
particles in a liquid film, when pi > 0.5. The precise form of m{pi) has a negligible effect 
on the qualitative behaviour of the system. Here we use 



a 1 
m{pi) = -{l + t^nh[30{pi --)]}. 



(17) 



Next we consider the time evolution of the hquid density pi(r, t). The dominant process 
governing the dynamics of the hquid is the evaporation and condensation of the hquid 
between the surface and the vapour reservoir above the substrate. We define two 
different chemical potentials: fi is the chemical potential of the liquid in the reservoir 
(cf. Eq. (4)) and iJ.s{r,t) = + /t* denotes the local chemical potential of the liquid on 
the substrate. We assume that the evaporative contribution to the time evolution of the 
liquid density is proportional to the difference between fis{^, t) and /i. This gives us the 
following expression: 



-Ml 



(18) 



dt ' 6pi{v,t) 

where the dynamical coefficient M"*^ is assumed to be a constant. The value of Mf^ 
determines the rate of the non-conserved (evaporation) dynamics of the liquid. We 
should also allow for (conserved) diffusive motion of the liquid over the surface. We 
assume from DDFT that the diffusion of the liquid takes a similar form as the diffusion 



of the nano-particles given in Eq. (14). We therefore model the full liquid dynamics by 



combining the diffusive and the evaporative terms 



dpi{Y,t) 
dt 



= V- 



MfpiV 



5F[pn.,pl] 

Spi{r,t) 



Mr 



.SF[pn,pi] 
Spi{r,t) 



(19) 
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The mobility coefficient Mi for the conserved part of the dynamics is assumed to be 
constant. The ratio between the conserved and non-conserved mobihty coefficients 
determines the inffuence that the diffusive/evaporative terms have on the overaU 
dynamics of the hquid (i.e. Mp^/M^ ^ 1 corresponds to the case when the hquid 
dynamics are strongly dominated by evaporation and Mp^/M^ <C 1 corresponds to the 
case when liquid diffusion plays an important role in the dynamics). Thus, equations 
([S]), ([9]), (14) and (19), taken together, define our model equations, which govern the 
dynamics of the system. 

Note that when the liquid density p/(r,t) is a constant, the theory reduces to the 
DDFT developed by Marconi and Tarazona [39j , that may be obtained by approximating 
the Fokker-Planck equation for a system of Brownian particles with overdamped 
stochastic equations of motion |[39l |40l [411 [42l |43] . Equations similar to (19) can also be 
derived in the context of hydrodynamics. The resulting mesoscopic hydrodynamic thin 
film equations contain different mobilities and local energies [48l HSj . The combination 
of the diffusive and the evaporative terms can also be seen as a combination of a 
conserved Cahn-Hillard-type dynamics with a non-conserved AUen-Cahn-type dynamics 



4. Phase Behaviour 

4.1. One- component fluid 

It is important to understand the equilibrium behaviour of the fluid in our system as 
this gives us some insight into how the system behaves when it is out of equilibrium. Of 
particular importance is to determine what phases we may observe and the stability of 
these phases. Since we are modelling the evaporative dewetting of the liquid we initially 
seek parameter values which lead to a high density liquid phase (liquid film) coexisting 
with a low density phase ('dry' substrate). Employing a linear stability analysis, we 
calculate the spinodal curve, i.e., the limit of linear stability for an infinitely extended 
system. The spinodal curve is defined as the locus of points where the curvature of 
the free energy is zero, = 0, which is equivalent to the isothermal compressibility 
being zero [45j . Note that in this section we set pi = p, for simplicity. We also calculate 
the binodal curve, i.e., the coexisting density values for a system in equilibrium, by 
equating the chemical potentials, temperature and pressure in each of the coexisting 
phases. The area outside the binodal curve is a stable region where we see no phase 
separation. Inside the binodal curve we have phase separation in the thermodynamic 
limit. However, the linear stability of the fiuid depends on whether the curvature of free 
energy is positive or negative. When we have positive curvature (outside the spinodal 
curve), the system at this state point rests within a local minimum of the free energy, 
i.e., it is linearly but not absolutely stable. There is a free energy barrier that must be 
traversed to cross into the (absolutely stable) equilibrium phase. This is known as the 
metastable region, where local fiuctuations in the density (if sufficient in size) create 
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Figure 2. Phase diagrams showing the binodal (sohd Hne) and spinodal (dotted 
line) for the one component pure fluid. In (a) we plot the phase diagram in the 
temperature T versus density p plane and in (b) we display the same phase diagram 
in the temperature versus chemical potential fi plane. 



nucleation points for the phase transition to occur. When the curvature of the free 
energy is negative (inside the spinodal curve) there is no free energy barrier. This is the 
unstable region where we have spontaneous phase separation, i.e. where fluctuations in 
the densities spontaneously grow. One may also say that the homogenous fluid layer is 
linearly unstable to harmonic perturbations with certain wave numbers. 

We first consider the phase behaviour of the reduced case where we have a single 
component fiuid with no nano-particles (i.e. = 0), as shown in Fig. [2| We set the 
external potential (/>^(r) = which results in a uniform fiuid density p^(r) = p = ^, 
where N is the number of 'particles' of liquid (i.e. filled lattice sites) and A is the area 
of the system. From Eq. ^ we find that the Helmholtz free energy per unit area of the 
uniform system is 

f = j = kBT[plnp+{l-p) ln(l - p)] - 2eip\ (20) 
We define the Helmholtz free energy per 'particle' as 



F 

a=- = kBT 



lnp+^^-^ln(l-p) 
P 



2eip. (21) 



From this, we may calculate other thermodynamic quantities: the pressure P and 
chemical potential //, which are given by the following relations [45j: 

^=p^f?y (22) 



dp 

/^=« + p(|^)- (23) 

We can calculate the spinodal for the one component fiuid from the free energy Eq. ([20|) 
and the definition of the spinodal ^ = 0, giving us the following equation: 

— = 4p(l-p). (24) 
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Equations (21) and ( |22[ ) give the following expression for the pressure in the system: 

P = -keT ln{l - p) - 2eip^. (25) 

In order to simplify the task of determining the phase diagram, we may use the symmetry 
of the Hamiltonian ([T]). The Hamiltonian remains unchanged under the exchange 
^ (1 — /^), i.e., it has a symmetry between 'holes' (nearly dry substrate) and 'drops' 
(liquid layer). This means that for the one component fluid, the phase diagram is 
symmetric around the line p = |, i.e. for a phase on the binodal with a density of p = pi 
the coexisting phase has a density of p2 = (1 — Pi)- Using Eq. (25) and this symmetry of 
the Hamiltonian, we obtain the following expression for the density along the binodal, 
by equating the pressure in the two phases (Pi = P2)'. 
ksT ^ 2(2p-l) 
ei ln[p/(l-p)] 

The maximum on this curve is at p = | and corresponds to the critical temperature 
ksT/ci = 1. Below the critical temperature there are two solutions; these are the 
coexisting densities. The binodal and spinodal curves for the pure liquid are plotted in 
Fig. [2](a). Fig. [2](b) shows the value of the chemical potential along these curves. 

The spinodal region can also be calculated from the dynamical equations (14) and 
(19) employing a linear stability analysis. This allows us to determine the typical 
length scales of the density fluctuations in the liquid fllm which might exist during 
the evaporation process. To perform the linear stability analysis we consider a liquid 
density which varies in space and time pi = p(r,t). The free energy for the single 
component fluid (p^ = 0) is given by: 



dr 



(27) 



where the subscript on the liquid interaction variable ci is dropped for simplicity. The 
steady state solutions of the liquid dynamical equation Eq. ( 19) represent the equilibrium 
density conflgurations. There are several steady states for this system, (e.g. a density 
proflle containing a free interface between two co-existing densities with p = Pcoex) but 
here we consider the simplest steady state: a flat homogeneous fllm with a density p = po 
which is deflned by: 

=0^ (28) 



6p 



PO 



We consider small amplitude perturbations 5p from po of the form p = po + 5p 



IS 



Po + 0e , where the amplitude <C 1 is a small positive constant, k 
the wave number and /3 is the rate of growth/decay with time (for positive/negative 
values) of the perturbation. We substitute this expression for p into the dynamical 
equation (19) and expand in powers of 5p. Then taking just the leading order terms 
allows us to derive a simple expression for /3 which can be solved analytically, 
exp 

Of 



A Taylor series expansion of the functional derivative of the free energy (27) yields: 
5F 
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(29) 



PO 



Substituting this approximation for the functional derivative (29) into the dynamical 



equation (19) we obtain 
I35p = McV 



Po{ik 



dl 
dp 



+ 



dp'^ 



5p + iek 5p) 



PO 



PO 



dp'^ 



PO 



Using the definition of po [Eq. (28)] gives |^|^^ 
terms 0{5p^) we arrive at the expression for the growth rate 



5p + ek^5p - /i) + 0{5p^). (30) 
and neglecting second order 



p 



/5 



^2 £ 



(31) 



When /5 is positive, small perturbations from the steady state po with the wave number 
k grow in amplitude over time. Conversely, if (3 is negative then small perturbations 
decay. Since Mc, po, and e are all positive quantities, (5 is always negative 

(i.e. the fiuid is stable for all wave numbers k) when the second derivative §^1^^ is 

positive. However when ^\p^ is negative, we find that (3 is positive for small values 
of k and negative when k is large - i.e. the fiuid is unstable against long wavelength 
(small wave number) fiuctuations in density. This corresponds to the thermodynamic 
definition of the spinodal as previously discussed. We may define a critical wave number 
/Cc, as the wave number at which /3(/Cc) = 0. When M^c 7^ and §^1^^ < 0, the critical 
wave number kc is given by: 



kc = ±\ 



Id^f 



e dp'^ 



(32) 



PO 



The real system we are modelling is very large (L 3> 27r//Cc), which means fiuctuations 
can occur on the full spectrum of wave numbers. The mode with wave number k = kj^ 
which has the largest positive value for /3 grow the fastest. The wave number k^ 
corresponds to a typical length scale 27r//c^ that are visible during the early stages of 
spinodal decomposition. However, at later stages (beyond the linear stage) the length 
scale of the modulations is likely to deviate from this value as the pattern coarsens. 
For the purely evaporative case when Mc = 0, the maximum value of /3{k) occurs at 
= = 0, which means there is no typical length scale in the early stages of the 
evaporation process. However, in the purely diffusive case M^c = (and §^1^^ < 0), 
/3{k = 0) = due to mass conservation, and the maximum value of /3{k) occurs 
at a non-zero value of fc^, so the typical length scale 271 /k^ is visible during the 
spinodal decomposition-evaporation. In all other cases the following expression is both 
a necessary and sufiicient condition for the existence of the typical length scale (i.e. 
fc^T^O): 



Mr 



< 



6 dp^ 



(33) 



PO 
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Figure 3. There are three possible forms for (a) shows an example of each 

form and (b) shows the regions in parameter space where these different types of p(k) 
curves occur, for the case when Mnc/Mc = 1 and ksT = 1. 

If a typical length scale does exist, then the corresponding wave number fc^ is given by: 

Figure [3] displays (a) typical f3{k) curves for a one-component fluid and (b) the 
values of 1/e and po where each type of curve is found (for the case when M^c/Mc = 1 
and ksT = 1). Three distinct cases are displayed: Case (i) the solid line (red online) 
in (a) and the striped area (red online) in (b) display the case when one observes the 
growth of density fluctuations with a typical length scale. Case (ii), the dashed line 
(blue online) in (a) and the hashed area (blue online) in (b) correspond to the situation 
when the fluid is unstable for density fluctuations corresponding to small wave numbers 
k but no typical length scale is observed, because the fastest growing mode is the k = 
fluctuation. Case (iii), the dotted line (green online) in (a), corresponding to the (green) 
solid area in (b) , displays the situation when the liquid is linearly stable against density 
fluctuations with all wave numbers k. 

4>2. Binary mixture 

To determine the binodal and spinodal curves for the binary mixture is more challenging 
than for the one component system. There are many more parameters deflning the 
model. In particular, the behaviour of the system strongly depends on the ratios between 
the different interaction strengths e^, and e^i- Ref. [53j gives a good overview of the 
equilibrium bulk fluid phase behaviour for binary fluid mixtures and the different classes 
of phase diagrams that may be observed. Here, we only describe the influence of the 
chemical potential of the nano-particles /i^ on the densities in the coexisting phases and 
on the critical point. 

For two phases of a binary mixture to coexist in thermodynamic equilibrium, there 
are four conditions that must be satisfled. We denote these two phases as (i) the 
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Figure 4. Binodal curves for the binary mixture for varying fixed values of the 
chemical potential of the nano-particles /i^, for the case when e^/e/ = 0.43 and 
^nl/^l = 0.57. 



low density phase (LDP) or 'dry' substrate and (ii) the high density phase (HDP) or 
substrate covered by a colloidal film: 

J.LDP ^ j.HDP^ (35^ 

Mf"" = /^r", (36) 
M^^^ = (37) 

pLDP ^ pHDP^ ^3g^ 

where T is the temperature, fii and /x^ are the chemical potentials of the liquid and 
nano-particles respectively, P is the pressure and the superscript denotes the phase. 
The first of these equations is trivial to satisfy in our model. We may then fix the 



chemical potential of the nano-particles, to some value ry, and then solve equations (36), 



(38), f^n^^ = T] and jJi^^^ = rj for the four density values: pf^^, Pn^^ ^ pf^^ p^^^ . 
The curves of the coexisting density values (binodals) for the parameters Cni/ci = 0.57 
and Cji = 0.43 are displayed in Fig. [4} The density values calculated for the two phases 
meet at a critical temperature to form binodal curves similar to the one found for the 
one component system Fig. [2^. However, we find that the reflection symmetry w.r.t. 
the liquid density pi = 1/2 seen for the one component fluid is broken. In particular, 
the critical point of the liquid is no longer at = 1/2. The liquid binodal reduces to 
the one for the pure liquid as /x^ ^ ±oc. The curves for Pn ^ —oo and /x^ ^ oc are 
identical within our model due to the symmetry of the Hamiltonian ([T]) . We observe that 
the density of the nano-particles ^ becomes very small as the chemical potential 
decreases p^ ^ — oc. Conversely, the density becomes very large ^ 1 as the chemical 
potential increases p^ ^ oc. The critical point on the nano-particle binodal curve shifts 
in a similar manner to that of the liquid binodal. 

We now consider the linear stability of the fluid for the two component case. There 
are many different possible steady states for the binary system including some with 
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interfaces between coexisting phases. However, here we hmit ourselves to investigating 
the hnear stabihty of the simplest steady state: where both components have a uniform 
constant density over the surface. The uniform density of the nano-particle film is 



denoted by = and the density of the liquid is given by pi 
determined from the condition: 
5F 



Pi 



Spi 



= 0. 



where p^ is 



(39) 



We consider small amplitude perturbations in the density of both components from this 
steady state. We assume these perturbations to take the form: 

p, = p° + (^ = p° + 0e^''V*, (40) 



Pn = Pn + XSp = Pn + X0e 



ik-r JSt 



(41) 



where the amplitude |(/)| <C 1 and the parameter x is the ratio between the perturbation 
in the densities of the two components. The sign of x determines whether any 
instabilities are in-phase (positive) or anti-phase (negative) between the two coupled 
density fields. From the magnitude of x we can determine whether an instability is 
driven by the liquid component (|x| <C 1), the nano-particles (|x| ^ 1) or stems from 
the interaction between the two components (|x| = 0(1)). Making a Taylor series 
expansion of the functional derivative of the free energy with respect to the liquid 
density keeping only terms linear in Sp and substituting this expression into the 



liquid dynamical equation (19) we obtain 



/3 = -(M^p°P + M^j( 



Pi tPu 



Pi ^Pn 



nl 



(42) 



dpidprj 

Turning our attention now to the nano-particles: in order to linearise Eq. ( |14[ ), we 
must first examine the nano-particle mobility function M^(p/,p^) = p^m(p/), given by 



Eqs. (16) and (17) in our model. Making a Taylor series expansion of the function m{pi) 
we obtain 

m 



a 



= jo + liSp + 0{5p'), 



(43) 



where 70 for small values of pi {pi < 0.45) and 70 1 for large values of pi 
{pi > 0.55). Substituting a Taylor series expansion of the functional derivative, together 
with Eq. ( [43| ), into the time evolution equation for the nano-particles Eq. (14) we find 



X/3 = -ap°7oP(x 



+ 



When pi is small, 70 



dpi 

0. In consequence, x 



(44) 



and the Eqs. (42) and (44) reduce 



to the one of the one-component fiuid (with a local free energy that depends also on 
p^). We now address the case when pi > 0.55, and we therefore assume 70 = 1. The 



expressions for the time dependency of the amplitudes of the two density fields (Eq. (42) 



and Eq. (44)) can be solved simultaneously to determine /3 and x as a function of the 
wave number k. This allows us to determine the stability of the fiuid for different 
values of the system parameters. A fact that simplifies the analysis is that these two 



Modelling the evaporation of thin films of colloidal suspensions using DDFT 



15 



equations may be written in matrix form (similar to the case of two coupled mesoscopic 
hydrodynamic equations for dewetting two-layer films [54j): 

/5l^) = M- G(^), (45) 



where, 



1 



M = 



G = 



-ap^r 

-{Mlp^ie + Ml)^ 

k^^n + fnn k^^nl + fnl 
k^^nl + fnl k'^^l + fll 



and we have used the shorthand notation fij = Lo^o, where i, 7 = n,/ and where 

opi Opj Pi Pyi 

n denotes the nano-particles and / denotes the hquid. We can determine (3{k) using the 
foUowing expression for the eigenvalues of a 2 x 2 matrix: 



,„).IMM^±^IM^_|M.G|. ,46) 

We define critical wave numbers k = kc for the density fiuctuations in the two-component 



fluid as the wave numbers at which one of the two solutions of Eq. (46) is equal to zero. 



Below [Eqs. (49) and (50)] we derive explicit expressions for kc which can be used to 
determine the conditions for the linear stability of the two component fluid (i.e. the fluid 
is stable when there is no solution for kc and the function f3{k) < for all wave numbers 
k). Since the matrix M is diagonal and all diagonal elements are non-zero the inverse 



M exists, allowing us to rewrite Eq. (45) as a generalised eigenvalue problem 



(G - M- V) j^^j = 0. (47) 

Setting /3 = in this equation in order to flnd the critical wave numbers fcc, we flnd that 
the determinant |G| = 0, implying that 

ktienei - ell) + kUcnfll + ejnn - ^Cnlfnl) + fnnfll " fnl = 0- (48) 

In the special case when e^e/ = e^^, the critical wave number is given by: 



/^^ , / fnl fnnfll ^^^^ 

^nfll + ^ifnn " "^^nlfnl 

However, more generally, when e^e^ 7^ e^^, we have 



fc-^; -'^^^ (50) 



where 



b — ^nfll + ^ifnn " "^^nlfnU 

C = fnnfll - fnl' (51) 
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Figure 5. (a) - (f) Dispersion relations for the case when e^e/ > e^^, enfii + ^ifnn — 
'^^nifni > (i.e. in Eq. (51) a > 0, 6 > 0) and ksT = 1. The graphs in the top 
row show f3{k) and the graphs in the middle row show x(^)- (a) (d) the 
parameters are: 0.648, = 0.3, 1.25, 0.5, e^^ 0.6, /i -3.35, 

1,M^ = and M^^ = 1. In (b) and (e) the parameters are: = 0.648, = 0.3, = 
1.25, en = 0.5, e^^ = 0.6, /i = -3.35, = 1,M^ = and M^^ = 9. In (c) and (f) the 
parameters are: p^ 0.901, pJl 0.3, 1.4, 0.6, e^^ 0.8, /i -3.8, = 
1, = and M^^ = 1. (g) shows how the spinodal line shifts for increasing p^, where 
— ^ni ~ 0.25, Cn = and /c^T = 1. (h) shows how the spinodal line shifts as the 
value of a = CnQ — e^^ increases, when p^ = 0.3, = ei and ksT = 1. 



Note that the locus c = is the spinodal curve [55l IS]. We categorise the linear 
behaviour of the system by the signs of a and b in Eq. ( [sT] ) as they have a profound 
impact on the shape of the dispersion curves. In Fig. [5] we display for the case when 
a > and 6 > all the different possible /3{k) curves, together with the corresponding 
x{k)- From Eq. (|46]) we see that there are two branches for which we denote 

/3^{k) and /3~{k). The corresponding x{k) curves are denoted and x~ ^ respectively. 
The branch (red solid lines in Fig. [5](a) - (c)) corresponds to the second term in 
Eq. (46) being positive and the /3~ branch (dashed lines) corresponds to the second 
term being negative. Due to mass conservation of the nano-particles, there is always 
one /3(/c) branch that is zero at = 0, and x = at = for the curve that corresponds 
to the /3{k) that is not zero at A; = 0. If fu > then /3^{k = 0) = and x~{k = 0) = 0. 
Alternatively, if fu < then /3-(fc = 0) = and x^{k = 0) = 0. 

Inspecting Fig. [5| we find that there are three possible forms for /5(/c), similar to 
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the one component fluid case (Fig. [3|). In Fig[5|(a) there is a maximum in /3+ at A; 7^ 0. 
This indicates that the fluid is unstable and that one wiU observe the growth of density 
fluctuations with a typical length scale 27r//c^ during the early stages of the spinodal 
process. Here fc^ refers to the wave number at the maximum of (5^ . In the second case 
[Fig. [5](b)] the fluid is linearly unstable for density fluctuations with small wave numbers 
k as in Fig. [5](a), but the maximum in f3^{k) occurs at = 0, i.e. there is no typical 
length scale visible in the density proflles. In the last case, shown in Fig. [5](c), the fluid 
is stable for all wave numbers k. We can use Eq. (50) for the critical wave number kc to 
determine the stability of the system. Since a > and 6 > we observe that there is at 
most one positive root of which only exists if c < (i.e. inside the spinodal). The 
spinodal curve when a > and 6 > is plotted in Fig. [5](g) and (h). In (g) we show how 
the spinodal line moves upwards in l/ci as the nano-particle density is increased. In 
(h) we show how the shape of the spinodal changes as the value of the parameter a is 
increased. We flnd that larger values of a make the spinodal curve flatter. 

In the reduced case when a = 0, one can determine from Eqs. (51) that b is 
always be positive. We flnd the same three types of dispersion relations for this 
reduced case as we do for the a > 0, 6 > case. For the case when a > and 
6 < we observe that the fluid is unstable for all possible combinations of parameter 
values. When we have a < (where at least one of the components of the mixture 
is more attracted to the other component than to itself), we flnd that b must be 
positive. In this regime, we observe that as ^ 00^ /3^{k) 00 and (3~{k) —00. 
This indicates that the density fluctuations with an inflnitesimally small typical length 
scale will grow fastest. This behaviour corresponds to a mixture that exhibits micro- 
phase separation. This behaviour is common in block copolymer systems, where 
chemical bonding prevents macroscopic demixing and instead demixing on the nano- 
scale is witnessed [57l EEJ • Micro-phase separation is also observed in certain colloidal 
suspensions [59l EQl EH El E] • 

An important consideration which must be made with the binary system is whether 
'liquid-particle' phase separation can occur as well as the 'liquid-gas' (low density - high 
density) phase separation that we have already discussed. The former corresponds to 
the coexistence of two phases, both having a high density. In terms of our system this 
would mean that we have coexisting phases with a high liquid density (i.e. pi ^ 0.6) in 
each phase. The coexisting values depend upon the temperature T, chemical potential 
p and the interaction energies e/, and Cni- It is known that the following condition 
must be satisfled for 'liquid-particle' demixing to occur [53] : 

em < (52) 

The existence of liquid-particle phase separation in addition to the gas-liquid phase 
separation implies that for certain parameter values we may have three co-existing 
phases. This situation has the potential to lead to dramatic consequences for the pattern 
formation in our dynamical system. We return to this issue in Sec. 5.4 below. 
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5. Nonlinear Dynamics 

We now go beyond the linear analysis presented above and discuss numerical results for 
the fully non-linear time evolution. 



5.1. Numerical setup 

Recall that the dynamics of our model is governed by the coupled partial differential 
equations ( 14) and (19), together with the free energy functional Eq. ([s]). We numerically 



solve these non-linear partial differential equations using a finite difference scheme on 
a square lattice with grid spacing Ax = 1. The time step size At varies between 
simulations, as the stability of our numerical scheme depends strongly on the values of 
the parameters in the model. Central difference approximations are made for the partial 
derivatives with respect to space and forward difference approximations are made for the 
partial derivatives with respect to time. The Laplacian terms (V^) are approximated 
using the eight-neighbour discretisation [64j : 

where, ^ p^^ denotes a sum over the nearest neighbour lattice sites and ^ p^^^ 
denotes a sum over the next nearest neighbour lattice sites. Alternative approximations 
may be used [65]. However, the choice of Laplacian approximation has little effect on 
the qualitative behaviour of the system. 

Our numerical results show that as the various parameters in the model are varied, 
several different patterns are formed. We begin in Sec. |5.2| by considering the effect 
of changing the chemical potential p of the vapour reservoir. We then focus on the 
fingering instability. In Sec. |5.3| we discuss the effect of varying the parameter a in the 
nano-particle mobility and also the role which the conserved part of the liquid dynamics 
plays in the overall dynamics of the system, by varying the mobility coefficient 
Finally, in Sec. |5.4| we discuss the influence of 'liquid-particle' demixing at the receding 
front and how this affects the fingering mechanism. 

For the simulation results shown in Sec. 15.21 and Sec. 15.31 we set the interaction 
energies to ei = 1.4, = 0.6 and Cni = 0.8. Using the linear stability analysis we have 
shown that it is possible to obtain stable phases with these parameter values [Fig. [5](c), 



a > and 6 > in Eq. (51 )\ 



5.2. Influence of the vapour chemical potential p 

In section [4T] we have discussed how changing the value of /x, the chemical potential of 
the vapour reservoir above the surface, affects the structures displayed by the system 
as the pure liquid evaporatively dewets from the surface. Recall that as p is decreased 
below its coexistence value /icoex, the dewetting mechanism is at first via the nucleation 
of holes and then when p is further decreased, via spinodal dewetting. This sequence 
as p is decreased is also observed when there are nano-particles dispersed in the liquid. 
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although as previously discussed, the phase boundaries are shifted and there is the 
possibility of other (liquid-particle) phase transitions. Choosing a binary mixture with 
the parameter values fc^T = 1.0, = 0.6, = 1.4 and Cni = 0.8 gives a phase with a 
high liquid density coexisting with a low density of the liquid. We also set M[ = and 
Mp^ = 1 to allow us to initially focus on just the evaporative non-conserved dynamics of 
the liquid. We set the initial density profiles corresponding to a (high density) uniform 
film of liquid with density p/(x,t = 0) = 1 — 10~^ mixed with nano-particles having an 
average density of = 0.3. In order to allow for the growth of density fiuct nations 
when the system is (linearly) unstable, we add a small amplitude random noise to 
the density profile of the nano-particles. Thus the initial nano-particle density profile is 
p^(x, t = 0) = p^^ + 2A(y — 0.5), where y is a random real number uniformly distributed 
between and 1 and A is the magnitude of the noise. Without these random density 
fiuctuations the density profiles would remain uniform under the evolution of the DDFT 
with a film of liquid remaining. Our boundary conditions are periodic in all directions. 

The c = spinodal curve (as calculated in the previous section. Fig. [5]) defines the 
limit of stability, i.e. inside this line the fiuid becomes linearly unstable. Thus, the fiuid 
is unstable when < —3.869, where /3 = l/fc^T. The speed of the process increases 
with decreasing values of /3fi. For very low values of the evaporation process is so fast 
that we do not see any pattern formation in the nano-particles - the liquid evaporates 
too quickly for the nano-particles to diflFuse. For the parameters ei = 1.4, = 0.6, 
Cni = 0.8, a = 0.5, = and M/^^ = 1 this occurs when f3fi ^ —4.2. Fig. [6] shows the 
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Figure 6. Density profiles displaying evaporation via spinodal decomposition. The 
top row shows the liquid density profiles and the bottom row shows the nano-particle 
profiles at times t/ti = 7 (left), t/ti =8 (centre) and t/ti =9 (right), where ti = 
The system parameter values are: ksT = 1, = 1.4, = 0.6, Cni = 0.8, M^^ — 0, 



1, a = 0.5, = -4.08 and A 0.005. 
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particular case when = —4.08. We see that the hquid behaves in a similar manner 
to that of a single component fluid by spontaneously de wetting everywhere. Now the 
evaporation is slow enough for the nano-particles to move into areas with a high density 
of liquid during this evaporative process which creates a flne network structure [as shown 
in Fig. [6](e) and (f)]. However, this diflFusion of the nano-particles is limited as it is still 
a much slower process than the evaporation of the liquid. We observe that towards 
the end of the process small heaps of nano-particles are formed, where the density is 
signiflcantly larger. This effect is enhanced by the attraction between the nano particles 
> 0). Increasing the value of /x, we move from the linearly unstable (spinodal) region 
into the metastable region of the phase diagram (c.f. Fig. [2]). Note however, that the 
actual values of on the binodal and spinodal now differ from those in Fig. [2] because 
of the inclusion of nano-particles, c.f. Fig. [4] 

On increasing the chemical potential into the range —3.869 < ^ —3.8, the 
liquid fllm becomes metastable but may still evaporate through the nucleation and 
growth of holes. Fig. [7| shows the case when = —3.86. The nucleation is caused 
by the random fluctuations in the density distributions (the initial density proflles 
are deflned in a similar manner to the previous case: pi{x^t = 0) = 1 — 10~^, 
Pn{x^t = 0) = + 0.4(y — 0.5). The amount of noise used and the free energy 
'barrier' for forming a hole determines the probability of a nucleation event occurring. 
There is a critical hole radius W which can be determined from the free-energy of the 
system. If a hole is smaller than this critical radius then it will shrink and the liquid 
density will return to its bulk high density value in this region. However, if the size 
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Figure 7. Density profiles displaying nucleation and growth of holes which leads to 
the development of a network pattern. The top row shows the liquid density profiles 
and the bottom row shows the nano-particle profiles at times t/ti = 30 (left), t/ti = 200 
(centre) and t/ti = 800 (right). The system parameter values are: ksT = 1, = 1.4, 
Cn 0.6, Cni 0.8, Mf 0, M^^^ = 1, a = 0.5, = -3.86 and A 0.2. 
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Figure 8. The solid lines display the critical hole radius in units of the lattice grid 
spacing a versus the chemical potential /i for different nano-particle densities pn- The 
lines start at the lowest value of ji for which the system is linearly stable. The dotted 
lines show the value of the chemical potential at coexistence /icoex, which the curves 
approach asymptotically. Here ksT = 1, = 1.4, = 0.6 and Cni = 0.8. 




of the hole is larger than this critical value then this hole will begin to grow. We can 
apply classical nucleation theory to calculate an estimate for the critical hole radius 
by determining the change in free energy AF when a low density (thin film) circular 
'hole' with a radius of R is inserted into the metastable liquid film. The radius R which 
corresponds to the maximum in AF is the critical hole radius R^. We approximate the 
change in free energy using the formula: 

AF = TiR^AP + 2^i?7, (54) 

where AF is the pressure difference between the two phases (the hole and the fiuid film). 
7 is the interfacial tension (excess free energy) for creating a straight interface between 
the two phases at coexistence /x = /Xcoex- It is important to note that the density values 
at coexistence are different to the density values out of coexistence. The density of the 
thin film of fiuid inside the hole is such that its chemical potential is equal to that of 
the bulk film of fiuid surrounding the hole. The critical hole radius R^ is given by the 
maximum of Eq. (54) - i.e. when = 0. Thus, the critical hole radius is F^ = 
Fig. |8] shows how R^ depends on the chemical potential ji for different nano-particle 
densities = 0, 0.1 and 0.3. This analysis only applies to the metastable region of the 
phase diagram (c.f. Fig. [2]) and therefore the critical hole radius R^ curves are bounded 
on the left by the spinodal curve and on the right by the binodal curve. Note that 
classical nucleation theory incorrectly predicts a finite value for F^ at the spinodal, due 
to the fact that the theory assumes a sharp interface as one approaches the spinodal. For 
further discussion on this see e.g. Ref. [66]. The probability for a hole to be nucleated by 
random thermal fiuctuations is proportional to e~^^-^ . The curves in Fig. [sjshow that 
as we approach /Xcoex the size of the critical hole increases. Hence, the probability of 
nucleation is greater nearer the limit of stability (spinodal curve) and decreases greatly 
as we approach coexistence (binodal curve). We observe that increasing the density of 
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the nano-particles Pn shifts the metastable region to lower chemical potential values and 
increases the range of the metastable region. This is due to the increase in the critical 
temperature associated with the increase in the nano-particle density p^, as previously 
discussed. 

In the case shown in Fig. [7| we are near the limit of stability where the critical 
radius of a hole is very small. This results in many nucleation points where holes are 
formed and begin to grow. The nano-particles are picked up by these growing holes 
which creates a rim around each hole with a high density of nano-particles in the rim. 
The holes in the liquid film continue to grow until their rims meet, creating a random 
polygonal network pattern of nano-particles. The liquid wets the surface of the nano- 
particles, which means the liquid remains on the surface in areas with a high density 
of nano-particles. This is due to the positive interaction energy between the liquid and 
the nano-particles {e^i > 0). 

If we increase the chemical potential further into the metastable range —3.8 ^ 
Pfi < /Xcoex, then the probability of a hole being nucleated becomes much smaller. In an 
experiment, in this parameter range, all holes that are formed are normally nucleated at 
defects or impurities of the surface (heterogeneous nucleation). Any interfaces between 
a high density liquid phase and a low liquid density phase will recede as the liquid 
evaporates. The velocity of the receding front depends on the value of /3iLi. For small 
values of we have a fast front, but the speed of the front reduces as /x approaches 
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Figure 9. Density profiles displaying the growth of an artificially nucleated point 
which develops branched structures. The top row shows the liquid density profiles 
and the bottom row shows the nano-particle profiles at times t/ti = 2000 (left), 
t/ti = 10000 (centre) and t/ti — 30000 (right). The system parameter values are: 



kBT = 1, = 1.4, en = 0.6, e^i = 0.8, Mf = 0, M^^ = 1, a = 0.5, 



-3.8 and 



A = 0.1. Note that instead of using Eq. ([53|, here the Laplacian term is approximated 
using: W = eihw + E P^^^ " 20p) . 
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/icoex- If we choose ji = /icoex then any straight front remains stationary. On a 
completely structureless substrate, one would usually expect such an interface to recede 
homogeneously; this is certainly the case for the pure liquid, when = 0. However, 
in our system when > we see the formation of fingers as the front recedes, due to 
the presence of the nano-particles. Fig. [9] shows a case when = —3.8. The initial 
density profiles in this situation differ slightly from the previous cases. Here we create 
an artificial nucleation point by setting the density of the liquid and the nano-particles 
to pi = Pn = 10~^ in a central 2a x 2a region. Without this seed nucleus the initial noise 
on the density profiles slowly decays and the densities of the two species return to their 
(metastable) equilibrium values. The liquid surrounding this nucleation point slowly 
recedes creating a circular dewetting front. As the front recedes, it begins to collect the 
nano-particles, as was also observed for the case in Fig. [7| However, here the growing 
hole does not meet any other holes and there is time for an instability to develop at the 
front which causes the liquid to evaporate faster in some regions and slower in others 
creating a 'wavy' front, as seen in Fig. |9](a) and Fig. [9](d). The 'bumps' at the front 
then appear to stop moving while the rest of the front continues to recede. As the front 
recedes and the hole circumference increases, more fingers develop, leaving a branched 
'fingered' nano-particle structure behind - Figs. [9](e) and (f). The time scale for this 
dewetting process is rather long and so we also observe some long-time coarsening effects 
on the finger structures. 

Recall that one of the goals of our work is to develop an understanding of how 
the different self-organised structures of nano-particles observed in the experiments 
j22l [251 [26l [23] are formed. Distinct observed structures are a) network structures and 
b) branched structures. Results from our model have shown how two different types 
of network structures can develop: i) a fine network structure created by a spinodal 
evaporation process (Fig. [g]) in which the nano-particle density varies over a fairly small 
range 0.27 ^ Pn ^ 0.45, ii) a large well defined network structure created by the 
nucleation and growth of holes in the liquid (Fig. [?]), in which the nano-particle density 
varies over a large range 0.05 ^ ^ 0.9. Our model also shows how instabilities at the 
evaporative dewetting front can create branched structures for certain parameter values 
(Fig. |9]). Note that there is also evidence of early stages of the fingering instability in 
the nucleation case shown in Fig. [7[ There one can observe that small bumps begin 
to develop in the edges of some of the larger holes. We now consider the formation of 
the branched structures in more detail. In particular, we investigate the dependence of 
these 'fingered' structures on the parameters of the model. 

5.3. Influence of mobilities on the fingering 

To make a detailed investigation of the branched finger structures it is important to 
maximise the distance a front can recede. This allows us to obtain better statistics 
which is important due to the fact that solving the DDFT in the fingering regime on 
a large grid can be time consuming. To achieve this objective we create a straight 
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Figure 10. Nano-particle density profiles for calculations with (a) a = 0.1, (b) a = 1, 
(c) a = 2 and (d) a = 3. In (e) we display a plot showing the dependence of the 
number of fingers on the parameter value a. The parameter values are: UbT = 1, 
ei = 1.4, en = 0.6, Cni = 0.8, Mf = 0, Mp^ = 1, = -3.8, Ax = 1 and A = 0.1. 



dewetting front along the bottom edge of the system (i.e. we set the two density values 
to Pi = Pn = 10~^ for the first five horizontal lines). We also set no- flux boundary 
conditions at the top and the bottom, to prevent a dewetting front forming at the 
top. The periodic boundary conditions on the left and the right side of the system 
domain remain. This set-up also allows for easier analysis of the branched structures 
since we begin with an initially straight front. We define a measure for the average 
number of fingers (/) to be the average number of branched structures per unit length 
in the final density profile, after the dewetting front has reached the top of the system. 
To calculate this quantity we implemented an algorithm which counts the number of 
transitions between a high density of nano-particles and a low density of nano-particles 
on each horizontal line of the system. We then determine the number of fingers on a 
given horizontal line by dividing this value by two. We set a minimum and maximum 
line for a given set of final density profiles and calculate the average number of fingers 
between these two lines. This value is then divided by the size of the system to give a 
value that is independent of the system size. We have investigated how this measure is 
affected by the different parameters of the system. 

We begin by discussing the effect that varying the parameter a has on (/) (recall 
that a determines the mobility of the nano-particles in the liquid film). We use the 
same parameter values as above, with /3p = —3.8 and a varying between 0.1 and 3. 
Fig. [To] shows final nano-particle density profiles for several values of a and also a plot 
of (/) versus a, which is calculated from the average of five runs. We see that the value 
of a has a significant influence on the average finger number. Increasing the mobility 
of the nano-particles results in fewer fingers being developed, which is in (qualitative) 
agreement with the experiments [23j and the KMC model results [29]. The mobility of 
the nano-particles directly influences the speed of the receding front, so evaporation is 
much slower when a is small. When a is very small {a ^ 0.002) we find that the two 
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density fields become practically decoupled, with the liquid evaporating at high speed 
leaving the nano-particles behind as a homogeneous film of the initial density . 



(a) (b) 




Figure 11. Nano-particle density profiles for (a) Mf = 0, (b) Mf = 5, (c) Mf = 10 
and (d) Mf = 15. In (e) we display a plot showing how the number of fingers (/) 
depends on the value of Mf. The parameter values are: /c^T = 1, = 1.4, = 0.6, 
Cni = 0.8, M^^^ = 1, a = 1, = -3.8 and A 0.1. 



The effect of liquid diffusion over the surface has also been investigated by varying 
the liquid conserved mobility M^. Using the same parameter values as above, and setting 
a = 1, we display in Fig. 11 final nano-particle density profiles for varying M^^ from 
to 15 and also the average finger number (/) versus averaged over five runs. We see 
that the diffusive mobility of the liquid M^^ does affect the average number of fingers but 
to a much smaller extent than the mobility of the nano-particles. The average finger 
number generally increases as is increased. 



5.4' Influence of liquid-particle demixing on the fingering 

We now discuss the effect of possible liquid-particle phase separation on the front 
instability. Such a phase separation may occur near the front even for nano-particle 
concentrations inside the liquid film that are far smaller than the binodal value for 
liquid-particle phase separation. This occurs because as a dewetting front recedes it 
collects nano-particles (as previously discussed) and therefore increases the value of Pn 
near the front. For certain parameter values we find that if Pn increases above a certain 
threshold value, then liquid-particle phase separation occurs in the front region. The 
liquid separates into two liquid phases, a mobile one poor in colloids and a less mobile 
one rich in colloids. To investigate the resulting effects we set the interaction energies 



to €/ = 1.7, Cni = 1 and vary from to 1.2. Eq. (52) indicates that we should observe 
liquid-particle phase separation for > 0.3. We set the average nano-particle density 
to be low, p^^ = 0.1, so that there is no liquid-particle phase separation in the bulk of 
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the fluid fllm. For the value of the chemical potential that we use ji = —4, the fluid is 
linearly stable for all values of and leads to a relatively fast de wetting front. 



(e) 0.02 
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Figure 12. Nano-particle density profiles for (a) = 0.3, (b) = 0.6, (c) = 0.8 
and (d) = 1.1. In (e) we display a plot showing how the mean number of fingers 
(/) depends on the parameter e^. When > 0.6 we begin to observe droplets being 
deposited along with the branched structures (as shown in (c)). The parameter values 



for these calculations are: ksT = 1, ei = 1.7, e„ 



1.0, M,^ 



0, Mp'^ 



1, a = 0.5, 



-4, = 0.1 and A = 0.025. 



In Fig. [T2](a)-(d) we display typical final nano-particle density profiles for various 
and the average finger number (/) versus e^, averaged over five runs. As is increased 
from Cn = 0, we initially see a linear increase in the average number of fingers (/). This 
reaches a peak at ^ 0.6, after which we begin to see the development of droplets. 
When 0.6 ^ ^ 0.9 we observe separation between regions of high density of nano- 
particles and low density of nano-particles occurring locally at the dewetting front. The 
areas with a lower density of nano-particles form very thin fingers which quickly rupture 
into a series of droplets, whereas the areas with a higher density form thicker fingers 
which are much more stable as shown in Fig. [T2|(c). The sections of the front with a 
lower density of nano-particles recede faster than the rest of the front. This results in 
a 'doublon' pattern which has been observed in many different systems, e.g. in the 
thin film directional solidification of a nonfaceted cubic crystal [67j. As we increase 
further, for > 1 we observe that the dynamics at the front remains similar, however, 
now all the finger structures are very thin and therefore quickly break up into droplets. 
One also notices that the tendency to form side branches decreases with increasing 
Cn whereas the orientation of the branches becomes increasingly perpendicular to the 
receding front. We also observe an increase in the density of the nano-particles within the 
fingers/droplets as IS increased. This is due to the increased attraction between the 
nano-particles. These results agree qualitatively with the KMC model [29]. The KMC 
results also show an initial increase in the number of fingers followed by a transition from 
fingers to droplets as the interaction energy between the nano-particles is increased. 
Fig. 17 of Ref. [29] displays a plot of (/) versus which shows a similar trend as the 
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DDFT results displayed in Fig. [T2|(e). However, as it is a discrete stochastic model, the 
details of the transition in the way the branching occurs are less discernible than in the 
present DDFT model. 



6. Concluding remarks 

We have presented a DDFT based model for the evaporative dewetting of an ultrathin 
film of a colloidal suspension. We have derived an expression for the free energy of the 
system using a mean-field approximation for a coarse-grained Hamiltonian model ([T]) 
for the system. We have also derived dynamical equations which describe the diffusive 
dynamics of the solvent and of the colloids as well as the evaporation of the solvent. 
We have considered the equilibrium phase behaviour of the pure solvent and of the 
two-component fiuid and identified parameter ranges where unstable, metastable and 
stable phases exist. We then solved the coupled dynamical equations numerically to 
investigate the different dynamical pathways of the phase transition and the resulting 
self-organised patterns of the nano-particles. 

The model successfully describes the various self-organised structures found in 
experiments and is in qualitative agreement with the discrete stochastic KMC model 
[29j . Our numerical results show how nano-particle network structures can form either 
from a spinodal processes (Fig. [6]) or through the nucleation and growth of holes (Fig. [t]). 
We have also observed how branched structures develop from a fingering instability of 
the receding dewetting front (Fig. [9]). 

The transverse front instability results from a build-up of the nanoparticles close 
to the front as the solvent evaporates, when diffusion is too slow to disperse them. 
This slows down the front and renders it unstable. As a result, density fiuctuations 
along the front grow into an evolving fingering pattern. This transverse front instability 
can be considered to be a self-optimisation process which maintains the mean front 
velocity constant [29j (see also the discussion of this in the context of a similar front 
instability occurring in the dewetting of non- volatile polymer films f68]). One may also 
say that the constant average front velocity is maintained by depositing some of the 
nano-particles onto the dry substrate creating the branched structures. Experimental 
observations show that the branched structures found in the ultra-thin film behind the 
mesoscopic dewetting front are initiated from random nucleation sites. The holes which 
are nucleated then grow, initially creating circular dewetting fronts. We subsequently 
observe that the fingering instabilities and the development of branched structures form 
on the circular interfaces. Fig. [9] shows how these circular branched patterns develop 
from a single nucleation point; our numerical results bare a striking resemblance to the 
experimental AFM images of this phenomena [23j . 

We have studied the branched structures in greater detail using a planar geometry, 
i.e., by creating initially straight dewetting fronts. We have considered how the different 
mobilities affect the fingering. The nano-particle mobility in liquid films has a significant 



effect on the average number of fingers in the branched structure (Fig. 10). The finger 
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number decreases rapidly with increasing mobility in agreement with earlier KMC results 
[29] . This behaviour can be attributed to the lower-mobility of the nano-particles 
that hinders re-distribution by diffusion and also reduces the speed of the dewetting 
front. For the system to attain a higher front speed it must deposit nano-particles onto 
the surface at a greater rate. Therefore, if the mobility of the nano-particles is low 
this leads to the creation of more fingers because in this case the average distance a 
nano-particle has to travel to reach a finger is smaller. Increasing the mobility for the 
conserved diflFusive dynamics of the liquid has the opposite effect on the average number 
of fingers (Fig. 11). The complex relationship between the diffusion of the liquid and 
the average finger number is not yet fully understood. Our hypothesis is that increasing 
the mobility of the liquid results in an effective increase in the speed of the dewetting 
dynamics of the liquid for fixed nano-particle mobility. Thus, the mobility of the nano- 
particles becomes lower in comparison. The increased finger number then results from 
the increased mobility contrast, in agreement with the general instability mechanism 
laid out above. 

The basic front instability as described above is a purely dynamic effect and does 
not depend on particle-liquid and particle-particle attractive interactions that favour 
demixing of the liquid and the nanoparticles. However, beside this regime (that we 
call the 'transport regime') we have investigated how interactions that favour demixing 
infiuence the instability (Fig. 12). In general, when increasing the interaction energy 
between the nano-particles one increases the tendency towards liquid-particle demixing. 
However, this has no practical effect as long as the nano-particle concentration is low, 
so that it is outside of the two-phase region. This is normally the case for our initial 
densities. However, in the course of the evaporative dewetting the density increases 
close to the receding front. Increasing the interaction energy between the nano-particles 
causes demixing to occur close to the front (but not in the bulk film). The demixing 
makes the fingering instability stronger (we call this the 'demixing regime'). At first, one 
finds a linear increase in the average finger number with increasing interaction energy. 
At higher values of e^, when the localised phase separation sets in, the fingers become 
straight with less side branches, before finally lines of drops are emitted directly at the 
front. In this regime, the mean number of fingers is determined by the dynamics and 
the energetics of the system. 

The results we have obtained with our DDFT model confirm that jamming of 
discrete particles (as already discussed in Ref. [37]) is not a necessary factor for the 
fingering instability to occur. Our model is a continuum model with a diffusion constant 
that is independent of the nanoparticle concentration. The present two-dimensional 
DDFT model has several advantages over the two-dimensional KMC model j28i I29j : In 
particular, the early instability stages are more easy to discern without the background 
noise of the KMC. Furthermore, the underlying free energy may be employed to 
analyse the equilibrium phase behaviour in detail, in a similar manner to Ref. [53j. 
Many standard tools for the analysis of partial differential equations can be applied 
to the coupled evolution equations, such as, e.g., the linear stability analysis of the 
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homogeneous films. In the future, one may perform a hnear stabihty analysis for the 
receding straight front and also investigate steady state solutions as has been done for 
evaporating films of pure liquids [49]. There are many details that would merit further 



investigations such as, for example, the doublon structure mentioned in Section [5^ and 
its relation to such structures formed in directional solidification [69j . 

The present DDFT model does not include the effect of surface forces, i.e., 
wettability effects (substrate-film interactions). Therefore, a non- volatile liquid film 
would not dewet the substrate. This implies that an important avenue for future 
improvement is to incorporate wettability effects into the model. This could be done 
by making a mean-field approximation to derive an expression for the free energy for 
a fully three-dimensional KMC model [TOl [Tj (after incorporating substrate-particle 
and substrate-liquid interactions). The resulting three-dimensional DDFT could then 
either be used directly or be averaged perpendicularly to the substrate employing e.g., a 
long-wave approximation. Another possible option consists of combining a mesoscopic 
hydrodynamic approach, e.g., a thin film evolution equation (see [HU [72[ [49]) with 
elements of DDFT. For a brief discussion of a similar approach see Ref. [73] . 

As a final remark, we recall that in the present work we have only considered 
dewetting from homogeneous substrates. However, it is straightforward to include 
surface heterogeneities in our model via the external potentials (/>i(r) in Eq. ([s]). As 
future work, it would be interesting to study the infiuence of surface patterning on the 
finger formation displayed by the present system. 
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